/*  sdfhypbl.cpp  by K.Tsuru  */
// function ID 3300 DRADIX
/*****************************************************************************
SDouble class
It provides the hyperbolic functions.
 cosh x for func == COSH_CALC, sinh x for func == SINH_CALC
 pfSeries = CoshSeries or SinhSeries(function pointer)
Both functions have nearly the same algorism, then bounded in one program for
an easy maintenance.When only one of them(Sinh or Cosh) is used, in order to
avoid to link another one you must pass the function pointer of the series
function for a small |x|.
******************************************************************************/
#ifndef SN_H
#include "sn.h"
#endif

SDouble Hyperbolic(const SDouble& x, int func, SDouble (*pfSeries)(const SDouble& x)){
    if(!x.Sign(3300)){
        if(func == COSH_CALC) return 1.0;
        return 0.0;
    }

    const double XBySeries = 8.0; //When |x| is less than this value it uses the series.
    const uint figBySeries = 16u; //When the figures of 'x' is less than this value it uses the series.
    const double xMax = (double)DFIGURES*M_LN10*(double)DRADIX_EXP_MAX;  // = 301759.17...
//check the argument, e^x < DRADIX^DRADIX_EXP_MAX
    double xD = fabs( doubleD(x, 0) );
    if(xD > xMax) x.SetError(x.OVERFLOW_ERR, "Hyperbolic", 3300);

//When 'x' has a short and small value, it is fast to use the series.
    uint fig = x.Last() - x.First() + 1u;
    int xe = x.NetRdxExp(), series = 0;
    // |x| is less than XBySeries.
    if( (fig <= figBySeries) && (xD < XBySeries) ) series = 1;
//When |x| << 1.0, the calculation "sinh x" by "Exp(x)" yields a cancelation.
    if(xe < 0) series = 1;              // |x| < 1/DRADIX
    if(series) return (*pfSeries)(x);   //The series of "sinh" and "cosh" rapidly converges.

// |x| >= 1/DRADIX
    SDouble r, dr;

    r = Exp(Dabs(x));       //r = e^|x|
    long ef=(long)x.MaxSize() - (long)r.NetRdxExp();
    if(ef<=0) dr.SetZero(); //(1/r) is much less than r.
    else {
        RealSize C;
        C.SetEffFig(uint(ef));  //It can decrease the effective figures.
        dr = DReciprocal(r);    //dr = 1/(e^|x|)
        C.SetEffFig(0);
    }
    if(func == COSH_CALC) r += dr;  //cosh x
    else {                          //sinh x
        r -= dr;
// A bug was fixed on March 13, 2001.
        if(x.Sign() <0) r.ChangeSign(); //This statement was lacked.
    }
    r = DsDiv(r, 2);    //r /=2;
    return r;
}